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Abstract. We analyze the interplay of synchronization and structure evolution in an evolving network of 
phase oscillators. An initially random network is adaptively rewired according to the dynamical coherence 
of the oscillators, in order to enhance their mutual synchronization. We show that the evolving network 
reaches a small-world structure. Its clustering coefficient attains a maximum for an intermediate intensity 
of the coupling between oscillators, where a rich diversity of synchronized oscillator groups is observed. In 
the stationary state, these synchronized groups are directly associated with network clusters. 

PACS. 05.45.Xt Synchronization; coupled oscillators - 05.65.+b Self-organized systems 



1 Introduction 

The emergence of coherent dynamics in ensembles of in- 
teracting active elements is one of the basic manifestations 
of self-organization in complex systems. A large number of 
natural phenomena can be ascribed to the organized ac- 
tion of individual agents, whose joint function gives rise to 
collective signals. Specifically, periodic dynamics in living 
organisms at many levels (molecular, cellular, social) have 
early been suggested to originate in the synchronization 
of many elementary oscillations PUT 

The most significant instances of synchronization in 
natural systems are found in the realm of life sciences. 
Synchronized dynamics is essential to the function of cel- 
lular tissues such as in the brain and in the heart 0], 
and is known to characterize certain forms of behaviour in 
social animals, from insects to humans 5,6,7 . The obser- 
vation that synchronization is ubiquitous in the function- 
ing of biological systems immediately rises the question 
on the evolutionary mechanisms that may have given ori- 
gin to this form of self-organization. Several models have 
been explored where dynamical parameters are modified 
in response to "selection pressure" in the form of learning 
algorithms, in such a way that the function of the system 
evolves towards a specified goal [HIE]- 

In particular, Gong and Van Leeuwen have consid- 
ered an ensemble of attractively coupled chaotic maps 
whose interaction pattern, described by means of a net- 
work, evolves in such a way that the mutual synchroniza- 
tion of individual motions is enhanced . The algorithm 
favours the interaction between elements whose internal 
states are similar. In the present paper, we consider a 
similar form of adaptive evolution in a network of coupled 
non- identical phase oscillators m . In this kind of system, 
the heterogeneity of the ensemble competes with coupling 



against the emergence of coherent dynamics ^2] • We show 
that the evolution of the interaction network makes it pos- 
sible to partially synchronize groups of oscillators, for cou- 
pling intensities well below the synchronization threshold 
of globally coupled ensembles. Correspondingly, the net- 
work evolves towards a structure with relatively high clus- 
tering, approaching a pattern with small- world properties 
-as also found to occur for networks of chaotic maps ^U] . 
In contrast with the latter system, however, in oscillator 
ensembles clustering turns out to be maximal for an inter- 
mediate value of the interaction strength. We characterize 
the emerging dynamical and structural properties of the 
ensemble as functions of the coupling intensity. 



2 The model 

Our model consists of an ensemble of N coupled phase 
oscillators, whose individual evolution is given by 



N 



(i) 



i = 1, . . . , N, where u>i is the natural frequency of oscilla- 
tor i and r is the coupling strength. The weights Wij define 
the adjacency matrix of the interaction network: W%j = 1 
if oscillator i interacts with oscillator j, and otherwise. 
The number of neighbours of oscillator i is Mj — Y^j Wij . 
Interactions are symmetric, so that Wij = Wji and the 
network is a non-directed graph. It is assumed that this 
network is not disconnected. As explained in the following, 
in our model the interaction network changes with time. 

During the evolution, the network is quenched over 
time intervals of length T. Along each one of these inter- 
vals we calculate the average oscillation frequency of each 
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oscillator, 



1 

T 



t+T 



(f>i(t')dt'. 



(2) 



It is well known that, if the coupling constant r is suffi- 
ciently large, two oscillators i and j whose natural frequen- 
cies are close enough will asymptotically oscillate with the 
same average frequency, fli — flj for T — ► oo. The collec- 
tive manifestation of this long-term correlation between 
the dynamics of oscillator pairs belongs, precisely, to the 
class of synchronization phenomena addressed to in the 
Introduction. 

In our model, after each interval of length T has elapsed, 
the following adaptive mechanism is applied to make the 
network structure evolve. An oscillator i is chosen at ran- 
dom, and the values = \Oi — fij\ are calculated for all 
j 7^ i. We detect the oscillator ji for which 6^ is min- 
imum amongst all the We also detect, amongst the 
neighbours of i, the oscillator j% for which <5y 2 is maximal. 
If ji is one of the neighbours of oscillator i, the network is 
not changed. Otherwise, the network link between i and 
_?2 is replaced by a link between i an j\ . After this update, 
a new interval of length T begins, and the process is suc- 
cessively repeated until some kind of stationary state is 
reached. 

To avoid the unnatural situation where the average 
frequencies are compared with too high (machine) preci- 
sion, we apply the above update mechanism only when 
the involved quantities 6ij 1 and 8ij 2 are larger than a cer- 
tain threshold e. In particular, if the maximal difference 
between the average frequencies of oscillator i and all its 
neighbours is below e, no changes are made. 

The mechanism by which the interaction network is 
rewired has been designed as to favour the connection be- 
tween oscillators with similar average frequencies. In other 
words, for a given coupling intensity, its main effect on the 
collective dynamics of the oscillator ensemble is to enhance 
the possibility of synchronization. This aspect is studied 
in Section 13.11 Now, what is the effect on the structure 
of the interaction network? Which topological features in 
the connection pattern emerge from such mechanism? To 
answer these questions we analyze statistical properties of 
the network connectivity, as described in Section [3.21 



3 Numerical results 

In the numerical calculations, we consider an oscillator 
ensemble of size N = 100. The natural frequencies Wi are 
chosen at random from a Gaussian distribution with zero 
mean and unitary variance, g(oj) = exp (— u> 2 / 2) / ^/2tt . 
Average frequencies are discerned with a threshold e = 
10~ 3 . Initially, the neighbours of each oscillator are cho- 
sen at random from the whole ensemble, establishing a 
link between each oscillator pair with probability p — 0.12. 
For N — 100, this implies that each oscillator has some 12 
neighbours on the average, which insures that essentially 
in all realizations of the initial condition the interaction 
network is not disconnected. As the same time, the connec- 
tions are rather sparse. The initial phases 4>i(0) are drawn 
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Fig. 1. Order parameter z versus time t, for four different 
values of the coupling strength r. From bottom to top r = 0.4, 
1.0, 1.5 and 2.0. 



from a uniform distribution in [0, 2tt). The evolution equa- 
tions are integrated using a standard Euler scheme, with 
time step At = 10 -2 . A rewiring of the network is at- 
tempted every 10 3 time steps, so that T = 10. 



3.1 Synchronization properties 

As a first step to characterize the interplay between the 
collective dynamics of the oscillators and the evolution of 
their interaction network, we study the synchronization 
properties of the ensemble. We begin our analysis consid- 
ering the synchronization order parameter |1 lj 



Z(t) 



1 

N 



N 



i<t>i(t) 



(3) 



The time average of Z(t), 



t+T 



Z{t')dt', 



(4) 



is calculated over intervals of length T, when the network 
is quenched, in the same way as the average oscillation 
frequencies, Eq. ©. The order parameter ranges from 
z ~ N^ 1 / 2 for unsynchronized motion, to z ~ 1 when 
the oscillators become fully synchronized. In Fig. ^ we 
show the time evolution of z for four different values of 
the coupling strength, r = 0.4, 1.0, 1.5 and 2.0, averaged 
over 100 initial conditions. In all cases, the order parame- 
ter reaches a stationary value which grows as the coupling 
strength increases. 

The analysis of z reveals the existence of partially syn- 
chronized states, even for very small values of the coupling 
strength. In order to obtain further information on these 
states we compare the distribution of the average frequen- 
cies fli with that of the natural frequencies u>i for repre- 
sentative single realizations. Figure [21 shows snapshots of 
the average frequency of each oscillator as a function of its 
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Fig. 2. Average frequency Oi versus natural frequency uji, 
for three different values of the coupling strength r. Each dot 
represents an oscillator in an ensemble of size N = 100. The 
three plots are snapshots at time t = 5 x 10 4 . Straight lines are 
the diagonals Qi = cut. 



natural frequency for three values of the coupling strength, 
r = 0.4, 1.0, and 2.0. For comparison, we recall that in a 
globally coupled ensemble -where Wij = 1 for all i and j, 
and Mi = N- the synchronization threshold with the same 
distribution of natural frequencies is placed at r c 1.6. 
The three snapshots are taken at time t — 5x 10 4 = 50NT, 
when z has practically reached its stationary level. At this 
time, on the average, each oscillator has been chosen 50 
times at the updates of the network structure. 

In the plots of J?j versus a;,-, synchronization is revealed 
by the presence of horizontal arrays of dots, correspond- 
ing to oscillators with different natural frequencies which, 
through the effect of coupling, have attained the same 
average frequency. Already for r — 0.4, we note the for- 
mation of many small groups of synchronized oscillators, 
all over the frequency distribution. A number of oscilla- 
tors, however, do not belong to any group. For r = 1.0 
several groups have mutually collapsed, and the resulting 
aggregates are larger. Finally, just a few groups remain for 
r = 2.0, containing essentially all the oscillators. 

In our system, thus, synchronized groups are already 
observed for coupling intensities well below the synchro- 
nization threshold r c of a globally coupled ensemble, quoted 
above. This is an indication that the evolution of the 
network structure succeeds at creating and maintaining 
connections between oscillators which are more likely to 
become synchronized, even for low coupling intensities. 
In a globally coupled ensemble -and, more generally, in 
a randomly connected ensemble with moderate to high 
connectivity- synchronization would first involve those os- 
cillators in the centre of the Gaussian distribution of natu- 
ral frequencies, where their number is larger jl 1U12| . This 
would typically give rise to a single synchronized aggre- 
gate around the mean natural frequency, surrounded by 
a "cloud" of non-synchronized oscillators. Increasing the 
coupling strength, the aggregate would grow in size at 
the expense of the "cloud". In our system, on the other 
hand, the appearance of small synchronized groups for low 
interaction strength reveals a non-trivial structure in the 
underlying network. Due to the presence of several groups, 
we may expect that the collective behaviour of the ensem- 
ble is dynamically richer than in the cases where a single 
aggregate forms. 
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Fig. 3. Time evolution of the clustering C for five different 
values of the coupling strength r. 

3.2 Network structure 

The synchronization properties discussed above -in par- 
ticular, the formation of synchronized groups for low cou- 
pling intensities- suggests that a reasonable choice to quan- 
titatively characterize the structure of the interaction net- 
work is the clustering coefficient C JSj • We recall that the 
clustering coefficient is a topological property of a net- 
work which measures the average number of neighbours 
of a given node which are in turn mutual neighbours. It 
is defined as the ratio of the total number of triangles to 
the total number of connected triples in the network. In 
a recent work |14) . it was found that both in random and 
in scale-free networks, an increase in the clustering coef- 
ficient favors the formation of oscillator subpopulations 
synchronized at different frequencies. 

In our model, as the network evolves, C changes with 
time. Figure |21 shows the time evolution of the clustering 
coefficient for several values of the coupling strength. Each 
curve is an average over 20 realizations. Even for r = 0, 
in the absence of interaction, C shows a small (but fast) 
growth from its initial level (Co ~ 0.058). In this case, 
the rewiring mechanism tends to cluster those oscillators 
whose natural frequencies are close to each other. For r > 
0, however, the effect is much stronger. The clustering 
coefficient attains a long-time asymptotic value Coo which 
depends on the coupling strength r, and can reach up to 
seven times the initial value. 

The asymptotic clustering coefficient as a function of 
the coupling strength is shown in Fig. Interestingly 
enough, it exhibits a maximum at r max « 1, where it 
reaches the value C sw 0.4. For r > r max , C M declines, 
and seems to approach the value observed for r = for 
sufficiently large coupling intensities. The presence of a 
maximum in for an intermediate value of r indicates 
that, as the result of the rewiring process, the resulting 
network acquires a more complex (less random) structure 
when interactions between oscillators are neither too weak 
nor too strong. When r is small, after a few rewiring events 
which connect oscillators with similar natural frequencies 
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Fig. 4. Stationary value of the clustering as a function of Fi S- 5 - Mean distance d as a function of time, for three differ- 
the coupling strength r (dots). The dashed curve is a spline ent values of the coupling strength r. 
approximation added for clarity. 



-thus slightly increasing the clustering coefficient- the net- 
work evolution ceases. Due to the modest effect of syn- 
chronization at such low interaction intensities, there is 
no much advantage, with respect to the adaptive rewiring 
mechanism, in further increasing C . On the other hand, 
for large coupling intensities, a high synchronization level 
is rapidly attained even for random interaction patterns. 
Again, consequently, rewiring has little effect on the net- 
work structure. The complexity of the resulting structure 
is therefore maximal for moderate interaction strengths, 
where the emergence of coherent synchronized oscillations 
is already significant, but such that the collective dynam- 
ics is not too much organized as to make the network 
evolution unnecessary. 

As a next step in the analysis of the network struc- 
ture, we consider the mean distance d between all oscilla- 
tor pairs. The distance between two oscillators is defined 
as the number of network links along the path of minimal 
length joining them. Figure [S] shows the time evolution 
of d for three different values of the coupling strength, 
r = 0.4, 1.0 and 2.0. The three curves are averages over 100 
different initial conditions. While, in the considered time 
span, d grows monotonically for r = 0.4, the evolution 
is manifestly non-monotonic for larger coupling strengths. 
At short times, d exhibits a fast growth, reaches a max- 
imum, and then decays slowly. As r increases, the maxi- 
mum is reached at shorter times, and the peak becomes 
sharper. Together with the observation that asymptotic 
clustering Coo is large, the small values of d for long times 
indicate that the network has evolved, from its initial ran- 
dom structure, to a small-world pattern. 

Why is it that the mean distance reaches a maximum 
for intermediate times? In order to suggest an answer to 
this question, in Fig. we present snapshots of the net- 
work structure, together with the corresponding plots of 
the average frequency J2j versus natural frequency u>i , for 
r = 2.0 at three different times. Networks were drawn 
using the Pajek software, optimized for display with the 
Kamada-Kawai algorithm 17 . Initially, the network is 




Fig. 6. Snapshots of the network structure and the corre- 
sponding average frequency J2; versus natural frequency uji for 
r — 2.0, at three different times, corresponding to the initial 
condition, the maximum in the mean distance d (see Fig. [SJ, 
and the stationary state. 



random and f2i — LUi for all i. For t = 10 , the network 
is clearly divided into clusters, which naturally leads to a 
larger value of the mean distance d. In fact, at this time, d 
reaches its maximum (see Fig.[SJ). The plot of average fre- 
quencies shows the presence of well-defined synchronized 
groups. Inspection of the state of individual oscillators re- 
veals that there is a direct correspondence between the 
groups and network clusters. The largest cluster is formed 
by mutually synchronized oscillators with fii — 0, while 
oscillators with frequencies far from zero form smaller clus- 
ters. When the system has reached its stationary state, at 
t = 5 x 10 4 , six groups of mutually synchronized oscil- 
lators can be clearly distinguished in the plot of average 
frequencies. The network structure shows that clusters are 
now better interconnected, showing a more compact over- 
all structure. This implies that the mean distance d must 
have decreased to a value close to that corresponding to 
the initial random network. 
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Fig. 7. Degree distribution P(k) for two values of the cou- 
pling strength r. The dashed curve corresponds to the Poisson 
distribution of the random initial network. 



structure, and the relation between synchronization and 
structure. 

We have begun our study by considering the synchro- 
nization order parameter z. The analysis of the time evo- 
lution of z shows that it reaches a stationary value that 
grows as the coupling strength between oscillators increases. 
This reveals the existence of partially synchronized states. 
In order to obtain more information on these states, we 
have inspected snapshots of the mean frequency versus 
the natural frequency, for different values of the coupling 
strength. This revealed the presence of groups of oscilla- 
tors which had become synchronized even for very small 
values of the coupling strength. 

The presence of several synchronized groups signals to 
a non-trivial structure in the underlying network. In or- 
der to analyze the network structure we consider first the 
clustering coefficient C. The time evolution of C shows 
that it reaches a stationary value Coo that strongly de- 
pends on the coupling strength r. In fact, Coo presents a 



Thus, the maximum at intermediate times in the mean 
distance must be ascribed to a transient effect in the struc- 
tural organization of the network, as groups of synchro- 
nized oscillators develop and segregate from each other. 
At those times, the network is divided into poorly in- 
terconnected clusters and the mean distance is relatively 
large. Later on, further collapse of synchronized groups 
increases their interconnection, and d drops. Note that, 
although counterintuitive at first sight, this behaviour is 
not incompatible with the monotonous time growth of the 
clustering coefficient. 

Additional characterization of the resulting network 
structure is provided by the degree distribution P(k), which 
gives the frequency of nodes with exactly k neighbours. 
We have analyzed the stationary distribution P(k) for 
various coupling strengths and, in Fig. we show re- 
sults for two representative values of r. For comparison, 
the Poisson distribution of the random initial network 
p(k) — exp(— A)A fe /fc!, with A = 12, is shown as a dashed 
curve. As r grows, the distribution deviates from the Pois- 
sonian shape. Its maximum flattens and shifts to the right, 
while the small-A; range becomes more populated. Even- 
tually, the distribution acquires a bimodal shape, with a 
secondary peak at small fc, as shown in Fig. for r = 2. 
Individual analysis of connectivities reveals that this peak 
represent the contribution of oscillators in small synchro- 
nized groups, whose frequencies are relatively far from 
zero. Large connectivities, on the other hand, belong to 
oscillators in the largest group. 



4 Conclusions 



the rewiring process, the network reaches a more complex 
structure when the interactions between the oscillators are 
neither too weak nor too strong. 

It is interesting to compare the behavior of the station- 
ary value of the clustering with the results obtained with a 
very similar model for chaotic maps JU|- In that case the 
stationary value of C was either zero for small coupling 
strength, or saturated close to 0.6 for larger values of the 
coupling strength. 

As a next step in the analysis of the network structure, 
we have considered the mean distance d. The observation 
that d is small for long times, together with the fact that 
the stationary clustering Coo is large, show that the net- 
work evolves from a random initial structure to a small- 
world network. Through an inspection of individual os- 
cillators we have established a direct connection between 
synchronized groups and clusters in the network struc- 
ture. This makes it possible to explain the non-monotonic 
transient behavior of the mean distance, and the bimodal 
shape of the degree distribution for large values of the 
coupling strength. 

It is worth recalling the synchronization phenomenon 
of flashing fireflies, as described by J. and E. Buck. They 
observed that fireflies brought to their hotel room in Bang- 
kok "...first flew about. ..then settled down in small groups" 
and finally "...the flashing within each group became mu- 
tually synchronous" [Sj. The similarity with the behavior 
observed in our model suggests that the mechanism of en- 
hancing interactions between dynamical elements whose 
internal states are similar plays a key role in the adaptive 
emergence of coherent dynamics. 



We have presented a model of adaptive evolution in a 
system formed by coupled non-identical oscillators. The 
attractive interaction between oscillators is described by 
means of a dynamic network, that evolves by favouring 
the interactions between elements whose average frequen- 
cies are similar. We have analyzed in detail the synchro- 
nization properties of the model, the underlying network 
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